##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone/GCClone_gene.all_record.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/DriverChoose/GeneLOH"
    type <- "IGC"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file
type <- opt$type

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info <- data.frame(fread(sample_info))
dat_gene <- data.frame(fread(gene_list))
dat_ccf <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
## 关注的病理亚型
if(type == "IGC"){
    dat_info <- subset( dat_info , Type == "IM + IGC" | Type == "IM + IGC + DGC" )
}else if(type == "DGC"){
    dat_info <- subset( dat_info , Type == "IM + DGC" | Type == "IM + IGC + DGC" )
}else if(type == "All"){
    #dat_info <- subset( dat_info , Type == "IM + IGC + DGC" )
    dat_info <- dat_info
}


dat_ccf <- subset( dat_ccf , ID %in% dat_info$ID )

###########################################################################################

## 判断基因多少比例出现在Trunk、Pre_Private、Inv_Private
result_plot <- c()
gene_list <- dat_gene$Gene_Symbol

for( geneN in gene_list ){

    tmp <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN )

    if(nrow(tmp) > 0 ){
        tmp$LOH <- ifelse( tmp$minor_cn == 0 , "LOH" , "Non-LOH" )
        tmp <- subset( tmp , Class != "IM" )

        res_tmp <- c()
        for( sample in unique(tmp$ID) ) {
            tmp_use <- subset( tmp , ID == sample )
            tmp_use <- unique(tmp_use[,c("ID" , "CLS" , "LOH")])

            ## 一个人若发生多个突变，算最早的
            if( length(which(tmp_use$CLS=="clonal [share]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [share]" )
            }else if( length(which(tmp_use$CLS=="clonal [early]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [early]" )
            }else if( length(which(tmp_use$CLS=="clonal [NA]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [NA]" )
            }else if( length(which(tmp_use$CLS=="clonal [late]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [late]" )
            }else if( length(which(tmp_use$CLS=="subclonal")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "subclonal" )
            }
            res_tmp <- rbind(res_tmp , tmp_use)
        }

        res_tmp <- data.frame(table(res_tmp$CLS , res_tmp$LOH))

        if(1!=1){
            ## Early和Share的合并定义为早
            early_share_loh <- sum(res_tmp[res_tmp$Var2=="LOH" & res_tmp$Var1 %in% c("clonal [early]" , "clonal [share]"),"Freq"])
            early_share_loh <- data.frame( Var1 = "Share/Early" , Var2 = "LOH" , Freq = early_share_loh )
            early_share_nonloh <- sum(res_tmp[res_tmp$Var2=="Non-LOH" & res_tmp$Var1 %in% c("clonal [early]" , "clonal [share]"),"Freq"])
            early_share_nonloh <- data.frame( Var1 = "Share/Early" , Var2 = "Non-LOH" , Freq = early_share_nonloh )
            early_share <- rbind( early_share_loh , early_share_nonloh )

            ## 其余的时间的合并定义为相对较晚
            other_loh <- sum(res_tmp[res_tmp$Var2=="LOH" & !(res_tmp$Var1 %in% c("clonal [early]" , "clonal [share]")),"Freq"])
            other_loh <- data.frame( Var1 = "Non-Share/Early" , Var2 = "LOH" , Freq = other_loh )
            other_nonloh <- sum(res_tmp[res_tmp$Var2=="Non-LOH" & !(res_tmp$Var1 %in% c("clonal [early]" , "clonal [share]")),"Freq"])
            other_nonloh <- data.frame( Var1 = "Non-Share/Early" , Var2 = "Non-LOH" , Freq = other_nonloh )
            other <- rbind( other_loh , other_nonloh )

            res_tmp <- rbind( res_tmp , early_share , other )
        }

        ## 其余的时间的合并定义为相对较晚
        other_loh <- sum(res_tmp[res_tmp$Var2=="LOH" & !(res_tmp$Var1 %in% c("clonal [share]")),"Freq"])
        other_loh <- data.frame( Var1 = "Non-Share" , Var2 = "LOH" , Freq = other_loh )
        other_nonloh <- sum(res_tmp[res_tmp$Var2=="Non-LOH" & !(res_tmp$Var1 %in% c("clonal [share]")),"Freq"])
        other_nonloh <- data.frame( Var1 = "Non-Share" , Var2 = "Non-LOH" , Freq = other_nonloh )
        other <- rbind( other_loh , other_nonloh ) 

        res_tmp <- rbind( res_tmp , other )
        res_tmp <- res_tmp %>%
        group_by( Var1 ) %>%
        summarize( Var2 = Var2 , Freq = Freq ,  Ratio = Freq/sum(Freq) )
        res_tmp$Hugo_Symbol <- geneN
        result_plot <- rbind(result_plot , res_tmp)
    }
}

out_name <- paste0(out_path , "/Driver.",type,".LOH.tsv")
write.table( result_plot , out_name , row.names = F , sep = "\t" , quote = F )


###########################################################################################
## 基因分布堆叠图
dat <- result_plot
dat$Ratio[is.na(dat$Ratio)] <- 0
dat$value_text <- paste0( round(dat$Ratio , 2) * 100 , "%")
dat <- subset( dat , Var1 %in%  c("Non-Share" , "clonal [share]") )
dat$Var1 <- ifelse( dat$Var1 == "clonal [share]" , "Share" , "Non-Share" )

## 计算p值

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

result <- c()
for( geneN in unique(dat$Hugo_Symbol) ){
    tmp <- subset( dat , Hugo_Symbol == geneN )

    a <- subset( tmp , Var1 == "Non-Share" & Var2 == "LOH" )$Freq
    b <- subset( tmp , Var1 == "Non-Share" & Var2 == "Non-LOH" )$Freq
    c <- subset( tmp , Var1 == "Share" & Var2 == "LOH" )$Freq
    d <- subset( tmp , Var1 == "Share" & Var2 == "Non-LOH" )$Freq

    if(length(a)==0){a=0}
    if(length(b)==0){b=0}
    if(length(c)==0){c=0}
    if(length(d)==0){d=0}

    p <- fisher.test( matrix( c(a,b,c,d) , ncol = 2 ) )$p.value

    if( p < 0.001 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
    }

    tmp$p_text <- ""
    tmp$p_text[1] <- p_text
    result <- rbind( result , tmp )
}

result$Var2 <- factor( result$Var2 , levels = c("LOH" , "Non-LOH") )
result$Var1 <- factor( result$Var1 , levels = c("Share" , "Non-Share") )

show_gene <- unique(subset( result , Var1=="Share" )$Hugo_Symbol)
result <- subset( result , Hugo_Symbol %in% show_gene  )
#gene_order <- c(
#    "TP53" , "ARID1A" , "CDH1" , "APC" , 
#    "ERBB2" , "PIK3CA" , "RNF43" , "MAP2K7" ,
#    "MTRR" , "MUC6" , "CFTR" , "BMP6" , "GAL3ST3")
#result$Hugo_Symbol <- factor(result$Hugo_Symbol , levels = gene_order )

plot <- ggplot( data = result , aes( x = Var1 , y = Ratio , fill = Var2 ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    facet_grid(.~Hugo_Symbol) +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".LOH.pdf")
ggsave( out_name , plot , width = 10 , height = 5 )


result <- subset( result , Var1 == "Share" )
plot <- ggplot( data = result , aes( x = Hugo_Symbol , y = Ratio , fill = Var2 ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    #geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 6,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".Share.LOH.pdf")
ggsave( out_name , plot , width = 6 , height = 2.5 )


###########################################################################################
## 计算整体上促进GC发生的基因其LOH
dat_use <- subset( dat_ccf , Variant_Classification %in% Variant_Type)

all_class <- c("All" , "IGC" , "DGC")
dat_res <- c()
for( use_class in all_class ){

    if(use_class!="All"){
        dat_driver <- subset( dat_use , Hugo_Symbol %in% gene_list & Class == use_class )
        dat_nondriver <- subset( dat_use , !(Hugo_Symbol %in% gene_list) & Class == use_class )
    }else{
        dat_driver <- subset( dat_use , Hugo_Symbol %in% gene_list & Class != "IM" )
        dat_nondriver <- subset( dat_use , !(Hugo_Symbol %in% gene_list) & Class != "IM")
    }
    driver_loh <- length(which(dat_driver$minor_cn==0))
    driver_nonloh <- length(which(dat_driver$minor_cn!=0))

    nondriver_loh <- length(which(dat_nondriver$minor_cn==0))
    nondriver_nonloh <- length(which(dat_nondriver$minor_cn!=0))

    p <- fisher.test( matrix(c(driver_loh , driver_nonloh , nondriver_loh , nondriver_nonloh) , ncol = 2) )$p.value
    if( p < 0.001 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
    }

    dat_plot <- data.frame( 
        type = use_class ,
        class = c("GCDriver" , "GCDriver" , "OtherGene" , "OtherGene" ) ,
        loh = c("LOH" , "Non-LOH" , "LOH" , "Non-LOH") ,
        mutnum = c(driver_loh , driver_nonloh , nondriver_loh , nondriver_nonloh) ,
        Ratio = c(
            driver_loh/(driver_loh+driver_nonloh) , 
            driver_nonloh/(driver_loh+driver_nonloh) , 
            nondriver_loh/(nondriver_nonloh+nondriver_loh) , nondriver_nonloh/(nondriver_nonloh+nondriver_loh)
            )
    )
    dat_plot$p <- ""
    dat_plot$p_text <- ""
    dat_plot$p[1] <- p
    dat_plot$p_text[1] <- p_text

    dat_res <- rbind( dat_plot , dat_res )
}
dat_res$value_text <- paste0( round(dat_res$Ratio , 2) * 100 , "%")

plot <- ggplot( data = dat_res , aes( x = class , y = Ratio , fill = loh ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    facet_grid(.~type) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 6,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.Driver_Other.LOH.pdf")
ggsave( out_name , plot , width = 6 , height = 4 )
